program main

!  Analisador e estimador de erros de discretização
!  com base nos estimadores de Richardson e GCI

	implicit none

	real*8, parameter :: great  = 1.0d+11
	real*8, parameter :: limite = 1.0d+10
	integer :: variavel ! número da variável de interesse nos arquivos de dados
	integer :: nm  ! número de malhas
	integer :: i   ! índice das malhas (i=1 --> malha mais grossa)
	real*8  :: r21 ! razão de h entre as malhas grossa(2) e fina(1)
	real*8  :: r32 ! razão de h entre as malhas supergrossa(3) e grossa(2)
	integer :: pL  ! ordem assintótica do erro verdadeiro
	!real*8  :: T_exato ! solução analítica exata
	real*8  :: Fs = 3.0d0	! fator de segurança do estimador GCI
	character*50,dimension(:),allocatable :: arq_numerico ! nomes dos arquivos das soluções numéricas
	real*8,dimension(:),allocatable :: h       ! comprimento de cada elemento da malha
	integer,dimension(:),allocatable :: nx     ! numero de volumes reais na direcao axial
	integer,dimension(:),allocatable :: ny     ! numero de volumes reais na direcao radial
	real*8,dimension(:),allocatable :: phi     ! solução numérica em cada malha
	real*8,dimension(:),allocatable :: phi_pU   ! Th extrapolado com as malhas 1 a 3 e pU_h
	!real*8,dimension(:),allocatable :: Eh      ! erro verdadeiro de Th
	!real*8,dimension(:),allocatable :: pE_h    ! ordem efetiva do erro verdadeiro de Th
	!real*8,dimension(:),allocatable :: pE_i    ! ordem efetiva do erro verdadeiro de phi_pU
	real*8,dimension(:),allocatable :: Uri		! erro estimado de Th com as malhas 1 a 3 e pU_h
	real*8,dimension(:),allocatable :: Ugci_pU! erro estimado de phi baseado em pU
	real*8,dimension(:),allocatable :: Ugci_p	! erro estimado de phi baseado em p
	real*8,dimension(:),allocatable :: pU		! ordem aparente do erro estimado de Th
	integer, parameter :: in   = 10 ! número do arquivo de entrada (dados)
	integer, parameter :: out  = 11 ! número do arquivo de saída principal
	integer, parameter :: out2 = 12 ! número do arquivo de saída secundário
	character*50  :: caso           ! nome do arquivo de saída
	character*90  :: title          ! título do gráfico
	character*102 :: head           ! título do gráfico + dia
	character*12  :: dia            ! data da simulação
	character*8   :: hora           ! horário da simulação
	character*90  :: note_caso      ! notepad + caso
	!character*50 :: arq_exato ! nome do arquivo da solução exata
	character*50 :: arq_dados ! nome do arquivo de dados principal
	integer :: dos                  ! acessa prompt dos
	integer :: j, k ! auxiliares

	call get_data

	open(out, file = "./output_files/" // trim(adjustl(caso)) // '.dat')

	call data_e_hora (dia, hora)

	head = trim(title)//" "//trim(dia)

	call write_data

	call calculate

	call present_results

	call get_graphs

	close(out)

	note_caso = 'gedit ./output_files/' // trim(adjustl(caso)) // '.dat'

	dos = system(note_caso)

contains
!===============================================================================
subroutine get_data

!	leitura do nome do arquivo principal de dados do caso a analisar

	arq_dados = 'input_file.dat'

	! leitura dos dados principais
	note_caso = 'gedit ./input_files/' // trim(adjustl(arq_dados))
!	dos = system( note_caso ) ! lista arquivo
	open(in,file = "./input_files/" // trim(adjustl(arq_dados)))

	read(in,*) variavel
	read(in,*) caso
!	read(in,*) com_solucao_analitica
!	read(in,*) arq_exato
	read(in,*) pL
	read(in,*) nm

	! alocação de memória
	allocate (h(nm), nx(nm), ny(nm), phi(nm), pU(nm), Uri(nm), Ugci_pU(nm))
	allocate (arq_numerico(nm), phi_pU(nm), Ugci_p(nm))

	! leitura dos nomes dos arquivos com resultados numéricos
	do i = 1, nm
		read(in,*) arq_numerico(i)
	end do

	read(in,*) title

	close(in)

	! leitura dos resultados numéricos nas nm malhas
	do i = 1, nm
	open(in,file = "./input_files/" // trim(adjustl(arq_numerico(i))))
		read(in,*) h(i)
		read(in,*) nx(i)
		read(in,*) ny(i)
		do k = 1, variavel
			read(in,*) j, phi(i)
		end do
	end do

	close(in)

end subroutine get_data
!===============================================================================
subroutine write_data

50 format(3x, 'titulo = ', a, //, &
   3x, 'Dia = ', a12, 5x, 'Hora = ',a8, //, &
   i4, 19x, ' = numero da variavel analisada', //, &
   a23, ' = nome do arquivo de dados principal', //,  &
   i4, 19x, ' = ordem assintotica do erro verdadeiro', //, &
   i4 19x, ' = numero de malhas', /)
51 format(t4,'h', t27, 'log10(h)', t40, 'r21', t53, &
   'arquivos com resultados numéricos:' )
52 format(1pe23.14, 1pe13.4, 16x, a)
53 format(1pe23.14, 2(1pe13.4), 3x, a)

	write(out,50) trim(adjustl(title)), dia, hora, variavel, arq_dados, pL, nm

	write(out,51)

	i = 1
	write(out,52) h(i), dlog10(h(i)), trim(adjustl(arq_numerico(i)))

	do i = 2, nm
		r21 = h(i-1) / h(i)
		write(out,53) h(i), dlog10(h(i)), r21, trim(adjustl(arq_numerico(i)))
	end do
	write(out,53)

end subroutine write_data
!===============================================================================
subroutine data_e_hora (dia, hora)

	integer*4 :: var(8)	! data e hora
	character*2 :: aux1, aux2
	character*4 :: aux3
	character*50 :: aux
	character*20 :: vardate, vartime, varzone	! data e hora
	character*12, intent(out) :: dia 			! data
	character*8, intent(out) :: hora 			! horário

	call date_and_time(vardate, vartime, varzone, var)

	write(aux,*) var(3)		! dia, uso de arquivo interno
	aux1 = trim(adjustl(aux))
	write(aux,*) var(2)		! mês
	aux2 = trim(adjustl(aux))
	write(aux,*) var(1)		! ano
	aux3 = trim(adjustl(aux))

	dia = '('//trim(aux1)//'/'//trim(aux2)//'/'//aux3//')'

	write(aux,*) var(5)		! hora
	aux1 = trim(adjustl(aux))
	write(aux,*) var(6)		! minuto
	aux2 = trim(adjustl(aux))
	write(aux,*) var(7)		! segundo
	aux3 = trim(adjustl(aux))

	hora = trim(aux1)//':'//trim(aux2)//':'//aux3

end subroutine data_e_hora
!===============================================================================
subroutine calculate

	real*8 :: psi	! taxa de convergência
	real*8 :: del	! auxiliar
	real*8 :: p		! the lowest order

! 	get aparent order pU
	pU = great
	do i = 3, nm
		del = phi(i) - phi(i-1)
		if (del /= 0.0d0) then
			psi = (phi(i-1) - phi(i-2))/(phi(i) - phi(i-1))
		else
			psi = 0.0d0
		end if
		if (psi > 0.0d0) then
			r21 = h(i-1)/h(i)
			pU(i) = dlog10(psi)/dlog10(r21)
		end if
	end do

!	phi_pU = phi extrapolado com as malhas 1 a 3 e pU
	do i = 3, nm
		r21 = h(i-1) / h(i)
		if (pU(i) <= 0.0d0 .or. pU(i) > limite ) then
			phi_pU(i) = great
		else
			phi_pU(i) = phi(i) + (phi(i) - phi(i-1))/((r21**pU(i)) - 1.0d0)
		end if
	end do

!	get richardson estimator Uri
!  Uri = erro estimado de phi com as malhas 1 a 3 e pU
	do i = 3, nm
		Uri(i) = phi_pU(i) - phi(i)
	end do

!	get GCI estimator
!	Ugci_pU = erro estimado de phi com as malhas 1 a 2 e pU
	do i = 3, nm
		r21 = h(i-1)/h(i)
		if (phi_pU(i) > limite) then
			Ugci_pU(i) = great
		else
			Ugci_pU(i) = Fs*dabs(phi(i) - phi(i-1))/(r21**pU(i) - 1.0d0)
		end if
	end do

!	get GCI estimator
!	Ugci_p = erro estimado de phi com as malhas 1 a 2 e p
	do i = 2, nm
		r21 = h(i-1)/h(i)
		p = dmin1(pL*1.0d0, pU(i))
		Ugci_p(i) = Fs*dabs(phi(i) - phi(i-1))/(r21**p - 1.0d0)
	end do

end subroutine calculate
!===============================================================================
subroutine present_results

	real*8 :: aux1, aux2, aux3, aux4, aux5 ! variáveis auxiliares

61 format (/, t3,'h', t17,'pU(phi)', t40,'pU(phi) - pL')
62 format(1pe12.4, 2x, 1pe22.14, 2x, 1pe22.14)
64 format(1pe12.4,2x,1pe22.14,2x,1pe22.14)
66 format ('# log(h)', t17,'Ordem aparente (pU)', t41, 'Ordem assintotica (pL)')
71 format (/, t3,'h', t15, 'phi', t37, 'phi_pU')
72 format(1pe12.4, 10(1pe22.14))
73 format(1pe12.4, 3(1pe22.14), 22x, 3(1pe22.14), 44x, 1pe22.14)
74 format(1pe12.4, 6(1pe22.14), 22x, 3(1pe22.14))
81 format (/,  t3,'h', t15,'Ugci(phi,p)', t37, 'Ugci(phi,pU)', t59, 'Uri(phi)')
82 format(1pe12.4, 8(1pe22.14))
83 format(1pe12.4, 2(1pe22.14), 22x, 3(1pe22.14), 44x, 1pe22.14)
84 format(1pe12.4, 5(1pe22.14), 22x, 2(1pe22.14))
86 format('# log(h)', t15, 'log|Ugci(p)|', t37, 'log|Ugci(pU)|', t59, &
   'log|Uri(pU)|')
96 format(/, '$ n_x   $&$ n_y   $&$ h', t33, ' $&$  \phi', t54, '$&$  p_U', &
   t74, '$&$ U_{ri}(\phi,p_U)', t94, '$&$ U_{GCI}(\phi,p_U)$\\')
97 format('$ 'i5, ' $&$ ', i5,' $&$ ', f10.8, ' $&$ ', 1pe15.8, ' $', &
   3('&   não se aplica   '), ' \\')
98 format('$ 'i5, ' $&$ ', i5,' $&$ ', f10.8, ' $&$ ', 1pe15.8, &
   3(' $&$  ---           '), ' $ \\')
99 format('$ 'i5, ' $&$ ', i5, ' $&$ ', f10.8, ' $&$ ', 1pe15.8, &
   ' $&$ ', f15.5, ' $&$ ', 1pe15.8, ' $&$ ', 1pe15.8' $ \\')

! 	write orders
	write(out,61)

	write(out,62) h(1)

	do i = 3, nm
		if (pU(i) < limite) then
			write(out,64) h(i), pU(i), pU(i) - pL
		else
			write(out,62) h(i)
		end if
	end do

	open(out2,file='./output_files/' // trim(adjustl(caso)) // '-ordens.dat')
	write(out2,66)
	do i = 3, nm
		if (pU(i) < limite) then
			write(out2,62) dlog10(h(i)), pU(i), 1.0*pL
		else
			write(out2,62)
		end if
	end do
	close(out2)

	write(out,71)

	write(out,72) h(1), phi(1)
	write(out,72) h(2), phi(2)

	do i = 3, nm
		if (phi_pU(i) < limite) write(out,74) h(i), phi(i), phi_pU(i)
		if (phi_pU(i) > limite) write(out,73) h(i), phi(i)
	end do

!	write uncertainties
	write(out,81)

	write(out,82) h(1)
	write(out,82) h(2), Ugci_p(2)

	do i = 3, nm
		if (Ugci_pU(i) < limite) then
			write(out,84) h(i), Ugci_p(i), Ugci_pU(i), Uri(i)
		else
			write(out,84) h(i), Ugci_p(i)
		end if
	end do

	open(out2,file='./output_files/' // trim(adjustl(caso)) //'-estimadores.dat')
	write(out2,86)
	do i = 3, nm
		if (Uri(i) /= 0.0d0 .and. Uri(i) < limite) then
			write(out2,82) dlog10(h(i)), dlog10(dabs(Ugci_p(i))), &
				dlog10(dabs(Ugci_pU(i))), dlog10(dabs(Uri(i)))
		else
			write(out2,82) dlog10(h(i)), dlog10(dabs(Ugci_p(i)))
		end if
	end do
	close(out2)

!	escreve dados para gerar tabela
	write(out,96)
	write(out,97) nx(1), ny(1), h(1), phi(1)
	write(out,97) nx(2), ny(2), h(2), phi(2)
	do i = 3, nm
		if (pU(i) < limite .and. pU(i) > 0) then
			write(out,99) nx(i), ny(i), h(i), phi(i), pU(i)/10, Uri(i), Ugci_pU(i)
		else
			write(out,98) nx(i), ny(i), h(i), phi(i)
		end if
	end do
   write(out,*)

end subroutine present_results
!===============================================================================
subroutine get_graphs

94 format("set grid", /, &
   "set time", /, &
   "set style data linespoints", /, &
   "set key right top", /, &
   "set yrange[-1:5]", /, &
   "set title '", a, "'", /, &
   "set xlabel 'log (h)'", /, &
   "set ylabel 'ordem aparente'", /, &
   "plot './output_files/", a, "-ordens.dat' using 1:2 title 'p_L',\", /, &
   "     './output_files/", a, "-ordens.dat' using 1:3 title 'p_U(phi)'", /, &
   "set term postscript eps enhanced color", /, &
   "set out './output_files/", a, "-ordens.eps'", /, &
   "replot")
95 format("set grid", /, &
   "set time", /, &
   "set style data linespoints", /, &
   "set key right bottom", /, &
   "set title '", a, "'", /, &
   "set xlabel 'log (h)'", /, &
   "set ylabel 'logaritmo do modulo das estimativas dos erros'", /, &
   "plot './output_files/", a, &
   "-estimadores.dat' using 1:2 title 'U_{GCI}(p)',\", /, &
   "     './output_files/", a, &
   "-estimadores.dat' using 1:3 title 'U_{GCI}(p_U)',\", /, &
   "     './output_files/", a, &
   "-estimadores.dat' using 1:4 title 'U_{ri}(p_U)'", /, &
   "set term postscript eps enhanced color", /, &
   "set out './output_files/", a, "-estimadores.eps'", /, &
   "replot")

	open(20, file = "./output_files/" // trim(adjustl(caso)) // "-ordens.gnu")
	write(20,94) trim(adjustl(head)), trim(adjustl(caso)), trim(adjustl(caso)), &
		trim(adjustl(caso))
	close(20)

	dos = system('gnuplot ./output_files/' // trim(adjustl(caso))//'-ordens.gnu')

	open(20, file = "./output_files/" // trim(adjustl(caso))//"-estimadores.gnu")
	write(20,95) trim(adjustl(head)), trim(adjustl(caso)), trim(adjustl(caso)), &
		trim(adjustl(caso)), trim(adjustl(caso))
	close(20)

	dos = system('gnuplot ./output_files/' // trim(adjustl(caso)) // &
		'-estimadores.gnu')

end subroutine get_graphs
!===============================================================================
end program main
